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A NONLINEAR PHYSICS-BASED OPTIMAL CONTROL METHOD FOR 
MAGNETOSTRICTIVE ACTUATORS 


RALPH C. SMITH* 

Abstract. This paper addresses the development of a nonlinear optimal control methodology for mag- 
netostrictive actuators. At moderate to high drive levels, the output from these actuators is highly nonlinear 
and contains significant magnetic and magnetomechanical hysteresis. These dynamics must be accommo- 
dated by models and control laws to utilize the full capabilities of the actuators. A characterization based 
upon ferromagnetic mean field theory provides a model which accurately quantifies both transient and steady 
state actuator dynamics under a variety of operating conditions. The control method consists of a linear 
perturbation feedback law used in combination with an optimal open loop nonlinear control. The nonlinear 
control incorporates the hysteresis and nonlinearities inherent to the transducer and can be computed offline. 
The feedback control is constructed through linearization of the perturbed system about the optimal system 
and is efficient for online implementation. As demonstrated through numerical examples, the combined 
hybrid control is robust and can be readily implemented in linear PDE-based structural models. 

Key words. Nonlinear optimal control, perturbation control, magnetostrictive actuators 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. This paper addresses the development of a model- based nonlinear optimal control 
method for magnetostrictive actuators in structural applications. Such actuators utilize the ‘giant’ mag- 
netostrictive effects provided by certain rare-earth compounds to produce significant strains in response to 
applied magnetic fields. As discussed in [7, 11], one core material which has proven very effective under 
a variety of operating conditions is Terfenol-D. Actuators utilizing this material can generate mechanical 
strains on the order of 500 fistrain in the linear range and up to 1000 fistrain in the nonlinear range. 
The materials are capable of generating forces in excess of 125 Ibf with specific values highly dependent 
upon transducer design. Furthermore, the material provides a broadband response ranging from DC up to 
20 KHz [17]. In combination, these qualities provide Terfenol-D transducers with significant capabilities as 
controllers and vibration absorbers in industrial and heavy structural applications. Such transducers have 
also been employed as sonar transducers and precision micropositioners. 

The full utilization of magnetostrictive transducers in all such applications requires quantification of 
the transducer dynamics in response to various inputs. Magnetostrictive materials such as Terfenol exhibit 
inherent magnetic hysteresis which is significant at moderate to high drive levels. Furthermore, numerous 
investigations have demonstrated the stress and temperature sensitivity of the materials along with the non- 
linear behavior of elastic properties such as the Young’s modulus [8, 28, 37]. Finally, the magnetomechanical 
relation between input currents and output strains is nonlinear and displays significant hysteresis at high 
drive levels. All such hysteresis and nonlinear effects must be incorporated in both the transducer models 
and control laws to utilize the full capabilities of the actuators at high drive levels. 
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(ICASE), M/S 403, NASA Langley Research Center, Hampton, VA 23681-0001. This research was also supported by the Air 
Force Office of Scientific Research under the grant AFOSR F49620-95- 1-0236. 
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The magnetization model we employ is based upon an extension of the ferromagnetic mean field model 
of Jiles and Atherton [23, 24, 25, 33] while magnetostriction and hence strains are incorporated through a 
quadratic domain rotation model [23] . As demonstrated through validation experiments in [9] , this combined 
model quantifies transducer dynamics for a large variety of prestresses and drive levels. The model also 
quantifies asymmetric minor loops which makes it appropriate for control design in structural applications 
which involve multiple frequencies and transient dynamics. 

We concentrate here on linear structural models which incorporate this nonlinear actuator model. Such 
linear models are common in applications characterized by large forces but small displacements and provide 
a natural regime for initial development of a nonlinear control method which incorporates the actuator hys- 
teresis and nonlinearities. A nonlinear open loop control is constructed first through the application of finite 
dimensional optimal control theory. This control adequately incorporates the hysteresis and nonlinearities 
inherent to the actuator but is not robust with regard to perturbations in operating conditions. Such robust- 
ness is provided by an additional feedback control constructed through linearization about the unperturbed 
optimal open loop control. The hybrid control containing the nonlinear open loop and linear perturbation 
closed loop components is highly robust, efficient to implement, and utilizes the flexibility and accuracy of 
the nonlinear actuator model to provide the capability for attenuating transient and broadband dynamics. 

To place this control method in perspective, we briefly summarize existing control techniques for magne- 
tostrictive actuators. For low drive level control applications, linear models and control methods have proven 
suitable for both bulk magnetostrictive actuators [5, 7, 31, 32] and magnetostrictive particle actuators [29]. 
In a similar vein, the effects of hysteresis were also neglected and a linear law employed in [38] when designing 
a high precision magnetostrictive micropositioner. Such methods break down at moderate to high drive levels 
due to inherent hysteresis and nonlinearities [13]. For example, hysteresis provides a phase lag effect which 
will destabilize a system if unaccommodated. One technique for extending the linear range of transducer 
dynamics is based on the assumption that the underlying system is linear with nonlinear output harmonics 
acting as a disturbance. As demonstrated by Hall and Flatau [17] and Hodges and Sewell [19], feedback 
techniques can then be employed to reduce disturbances and improve linearity for certain operating regimes. 
Various nonlinear control techniques have also been employed for high drive level applications. Jenner et al 
[22] developed an active vibration controller for predescribed wave forms by considering a nonlinear control 
technique implemented through switching between positive and negative gains to the actuator. A similar 
objective was attained via neural network controllers by Bryant et al [4]. 

Control laws based upon Preisach models have also been employed for a variety of smart materials 
including magnetostrictives [14]. Such models are based upon polynomial or piecewise constant approxima- 
tions to the nonlinearities and hysteresis loop, and are advantageous when the underlying physics is not well 
understood or quantified. Such characterizations provide a control input operator which is easily inverted (or 
has an inverse which is easily approximated) which facilitates control design based upon output linearization 
[36]. Feedforward control methods based upon Preisach models followed by linearization have been employed 
for piezoceramics [15] and are applicable for magnetostrictives in certain regimes. 

The generality of Preisach models, which provides their advantage when the physics is not well un- 
derstood, also leads to inherent limitations in many control applications. Because such models are not 
physics- based, they typically do not provide the capability for adapting to changes in operating conditions 
(e.g., drive levels, prestresses, minor loops) through the monitoring of system inputs. The transducer dynam- 
ics must be known a priori and incorporated directly in Preisach models whereas physics-based models of the 
type employed here can adapt to changing dynamic levels solely through the measurement or designation of 
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input currents. This limits Preisach- based control laws to predefined trajectories (e.g., periodic) and does not 
provide the capability for directly attenuating unmodeled disturbances or inputs. Moreover, for unanticipated 
initial conditions, such methods also lack the capability for controlling transient dynamics. Finally, Preisach 
models typically require a large number of nonphysical parameters which limits their flexibility and increases 
implementation time in many applications. 

The development of the physics-based control method is presented as follows. Section 2 contains a 
brief description of a typical magnetostrictive transducer along with an outline of the energy- based model 
employed in [9]. The modeling and approximation of transducer inputs to a thin beam are presented in 
Section 3. This provides the prototypical control system with nonlinear actuator inputs. The control prob- 
lem is discussed in Section 4. Following an outline of nonlinear optimal control theory for finite dimensional 
systems, a linear optimal control method is considered. Numerical results demonstrate the success of the 
method at low drive levels and its failure at high drive levels due to unincorporated phase lag effects. An 
open loop nonlinear control method which fully incorporates material nonlinearities and hysteresis is then 
developed. Numerical examples are used to show that this method provides excellent attenuation when 
the system is known exactly but is not robust with respect to system uncertainties. The final subsection 
of Section 4 illustrates the development and performance of a perturbation feedback control method ob- 
tained through linearization about the optimal control system. This feedback method provides excellent 
attenuation of structural dynamics, is highly robust with respect to operating uncertainties and is feasible 
for implementation. Finally, the method is effective for systems exhibiting broadband responses and both 
periodic and transient dynamics. 

2. Magnetostrictive Actuator Model. The issues which must be addressed when developing a 
nonlinear modeling and control methodology are illustrated through consideration of the transducer depicted 
in Figure 1. This construction is typical for actuators currently employed in structural applications (see 
[16]), and its dynamics exhibit the full range of nonlinearities and hysteresis which must be characterized 
and incorporated in control design. 

The primary components of the transducer consist of a magnetostrictive Terfenol-D rod, a surrounding 
wound wire solenoid, a surrounding permanent magnet, and a prestress mechanism consisting of spring 
washers and/or compression bolts. The input to the actuator consists of a time- dependent current X(t) to 
the solenoid. This generates a magnetic field H and corresponding magnetic flux B and magnetization M 
within the Terfenol rod. The rod is constructed so as to contain a large number of regions in which moments 
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Figure 1 . Cross section of a typical Terfenol-D magnetostrictive transducer. 
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are aligned perpendicular to the longitudinal rod axis (the orientation of regions, termed domains, is further 
aligned by the prestress mechanism). The application of the magnetic field causes the rotation of these 
moments which in turn generates strains and forces within the material. This provides the mechanism for 
actuation. We note that the magnetostrictive materials can also be used for sensing through measurement 
of the magnetic fields generated by stress- induced domain rotations. 

As illustrated in Figure 2, the relationship between the applied field H and the induced magnetization M 
displays significant hysteresis and saturation effects at high drive levels. This implies that the permeability 
fi, which relates the two, is a nonlinear, multivalued map. The magnetomechanical effects are also nonlinear 
as illustrated in Figure 3. At moderate drive levels, the relationship between the magnetization M and strain 
e is approximately quadratic, as depicted in Figure 3a, which yields the ‘butterfly 1 relationship shown in 
Figure 3b (the asymmetry is due to the use of experimental field input data when computing the modeled 
strain). As illustrated in [9], the magnetization /strain relation also exhibits hysteresis at high drive levels 
which must be incorporated in the magnetomechanical model. 



Magnetization M x 1Q 5 Magnetic Field H % 1() 4 


(a) (b) 

Figure 3. (a) Input/output relationships; (a) magnetization M and strain e and (b) Magnetic field H and 
strain e. 
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The model described in [9] is used to characterize the transducer dynamics. The magnetization compo- 
nent of the model is based upon the Jiles- Atherton mean field theory for ferromagnetic materials [23, 24, 
25, 33], This theory is based on the quantification of energy losses due to domain wall intersections with 
inclusions or pinning sites within the material (the transition regions between domains are termed domain 
walls). For a material which is free from inclusions, the domain wall movement is reversible which leads to 
anhysteretic (hysteresis free) behavior. However, materials such as Terfenol contain second phase materials 
which impede domain wall movement. At low field levels, domain wall movement about pinning sites is 
reversible and yields a reversible magnetization M rev . At higher drive levels, domain walls intersect remote 
pinning sites which provides an irreversible component Mj rr . It is this latter component which incorporates 
the energy loss and hysteresis in the model. 

To characterize the total magnetization M, we consider first the effective field within the material. For 
rods subjected to a constant prestress <7o, the effective field is given by 

H eff (t) = H(t)+aM{t) 


where 


H(t) = nl(t) 


denotes the magnetic field generated by a solenoid having n turns per unit length with an input current 2(t). 
The parameter a quantifies magnetic and stress interactions. Through thermodynamic considerations, the 
anhysteretic magnetization is then defined in terms of the Langevin function 


( 2 . 1 ) 


M an (t) = M s 


coth 


( a 

V a ) \H ef f(t) 


Here M s denotes the saturation magnetization of the material and a is a parameter which characterizes the 
shape of the anhysteretic curve. Energy balancing (see [24]) is then used to quantify the irreversible and 
reversible magnetizations through the expressions 


dM irr _ dJ M an (t) - M irr (t) 
dt dt kS — o\M an (t) — Mi rr (t)\ 


and 


(2.3) M rev (t) = c[M an (t) - M irr (t)} 

(J = ±1 while the constants c and k are estimated from the experimental hysteresis curves). Finally, the 
total magnetization is given by 

(2.4) M(t) = M rev (t) + M irr (t ) . 

To first approximation, the strains generated by the material are given by the bulk magnetostriction 

(2.5) A(t) = ^M 2 (t) 

where X 3 denotes the saturation magnetostriction (see [23] for details). In combination, (2.1)~(2.5) charac- 
terize the relationship between the input current 1 and the strains generated by the transducer. Details 
regarding the well-posedness of the model are given in [35]. 
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3. Structural Model. A structural system which has been experimentally employed to ascertain ca- 
pabilities and properties of magnetostrictive transducers (see [10]) is illustrated in Figure 4. This system 
consisted of a cantilever beam with end- mo unted actuators. Diametrically out-of- phase currents to the actu- 
ators generated bending moments which were used to attenuate transverse beam vibrations. We will employ 
a model for this experimental setup as a prototype to illustrate the optimal control method proposed here. 



Figure 4. Cantilever beam with magnetostrictive actuators. Uniform force inputs are depicted above the 
beam while the measurement point is indicated by the lower arrow. 


3.1. Thin Beam Model with Nonlinear Actuators. For modeling purposes, the beam is assumed 
to have length E, width b, and thickness h. The density, Young’s modulus, Kelvin- Voigt damping coefficient 
and air damping coefficient for the beam are denoted by Eb, co h and 7, respectively. The cross-sectional 
area of the Terfenol rod is denoted by A ma5 while the Young’s modulus and damping coefficient for the 
Terfcnol rod are denoted by E H and . The length and width of the connecting bar are denoted by E r and 
b r , respectively, while the bar density is given by p r . Finally, the transverse beam displacement is given by 
w while g(t, x) denotes an exogenous surface force to the beam. 

Moment and force balancing yields the strong form of the Euler-Bernoulli equations 


d 2 w dw . 

PW-Qjz&x) + T'Tft (*’ X ) + 


dx 2 


(t,x) = g(t,x) + 


d 2 M mag 

dx 2 


(t,x) 


0 < x < t 
t > 0 


™(*,0) = ^-(t,0) = 0 


t> 0 , 


along with appropriate initial conditions, as a model for characterizing the transverse beam dynamics. As 
detailed in [34], the composite density and internal bending moment are given by 


p(*E) — Pbhb + 2prb r ErXrod{x) 

M int (t,x) = EI(x)^^(t,x) + c D I-^^(t,x) 
where the characteristic function Xrod delineates the location of the rods and 


EI{x) = + 2A mag E H (h/2 + 4) 2 Xrod(x) 

c D I(x) = Cp ^ - + 2 A mag CD ( h/2 + 4) 2 Xr od(x) . 
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For the case when the Terfenol rods are driven diametrically out-of-phase, the external moment is derived 
from (2.5) and is given by 


M mag {t,x) = + 2M(t)M s } X rod{x) 

where /C M = (3X s /M 2 )A mag E H (h/2 + i r ) 2 - The inclusion of the weighted magnetization 2 M(t)M s provides 
the bias necessary to attain bidirectional strains. 

To obtain a weak form of the model, we take the state to be the displacement w in the state space 
X = L 2 ( 0,^) with the inner product 



p(j)tp dx . 


The space of test functions is taken to be V = Hl( 0,£) = {<f> e H 2 (0,£) | 0(0) = 0'(O) = 0} with the inner 
product 


<0,0)v = [* Eltff'l/'dx. 

Jo 

It should be noted that with these choices, V is continuously and densely embedded in H . Hence one has 
the Gelfand triple 


F<— X~X* ^V* 


with the pivot space X. 

A weak form of the model is then given by 


(3.1) 


I pivipdx+ I r yw r ipdx+ I Mint^” dx 
Jo Jo Jo 



Xi.fjiag'lp dx “h 



for all 0 £ V. It is in this form that we develop the approximation method and formulate the control 
problem. 


3.2. Approximation Method. A necessary step for constructing an implementable control law is the 
approximation of the infinite dimensional system (3.1). We employ a Galerkin approximation in the spatial 
variable to obtain a semidiscrete ODE system in time which is amenable to control formulation. Specifically, 
the spatial basis is taken to be where denotes the j th cubic B-spline modified to satisfy the 

fixed left boundary condition. Approximate solutions 

m+ 1 

(3.2) w m (t,x) = ^ Wj{t)4>j{x) 

3 = i 

are then considered in the subpace V m = span {(f) j}. To obtain a vector ODE system, the infinite dimensional 
system (3.1) is restricted to V m and posed in first-order form to yield 


(3.3) 


y(t) = Ay(t) + [B(u)](f) + G(t) 
2/(0) = 2/o • 
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The component system matrices have the form 


A = 


0 I 
Q~ l K Q~ l C 


[£(«)](<) = [M 2 (u) 4 -2M(u)M,\ (t) 


0 

Q-'B 


G(t) = 


Q-'m 


where y{t) = [u>i(f), • • • , • • • , it> m +i(f)] and 


(3.4) 


[Q)ij = f pMj dx [B]i = K m [ <j>" dx 

Jo J mag 

[K) t] = f Eltft'j dx [p(t)]i = f g(t, x)<t>i dx 
Jo Jo 

\C\ij = [ c d I4)"4>" dx . 

Jo 


Note that u(t) — X(t) denotes the control input to the system. For future development, it is useful to let b 
denote the 2 (m 4- 1) vector b = [0 , Q~ l B) T so that the control input can be written as 


(3.5) 


[B(«)](t) = [M 2 (u) + 2 M(u)M.](t) b. 


The system (3.3) provides the constraints employed in the control problem. 

3.3* System Parameters. For the examples which follow, the choice m = 12 was sufficient for resolv- 
ing beam dynamics in the frequency range considered and all reported results were obtained with m — 16. 
The dimension of the state vector y was then 34 x 1 due to the inclusion of both displacement and velocity 
components. 

The specific physical parameters employed in the examples are summarized in Table 1. It should be 
noted that the beam parameters are consistent with typical values for al uminum laboratory beams while the 
Terfenol parameters are within the range obtained for model fits to an experimental transducer [9] . For this 
choice of beam parameters, the first two natural frequencies for the system occur at 6.1 Hz and 38.3 Hz. 
To account for the effects of parameter discontinuities due to the actuators and damping in the system, it 
was necessary to obtain these values through a fast Fourier transform (FFT) of time domain data resulting 
from a simulated impact to the beam (it is not possible to obtain analytic expressions through separation 
of variables). The driving frequency in the numerical examples will be chosen close to but not exactly 
concurrent with these natural frequencies. 
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Beam 


Actuator 


Terfenol 


E b = 7.0861 x 10 10 N/m 2 
pb — 2863 kg/m 3 
c Db = 9.3663 x 10 5 Ns/m 2 
7 = .013 Ns/m 2 


E h = 7.0 x 10 10 N/m 
p r ~ 8524 kg/m 3 
i r = -0254 m 
b r ~ .002 m 
Amag — .0064 m 2 


a = 7105 A/m 
k = 7002 A/m 
a = .007781 
c= 0.3931 

M s = 1.3236 x 10 5 A/m 
X 9 — 9.96 x 10“ 4 


Table 1. Parameters for the beam and Terfenol transducer. 


4. Control Problem. The general form of the finite dimensional control system under consideration 


is 


(4.1) 


y{t) = f(y(t),u(t),t) 
y(t o) = yo 

with states y(t) G IR 2(m ’ 1 1 and controls u(t) € IR P where p - 1 for the case of a single actuator pair. As 
detailed in [6, 26, 27, 30], an appropriate performance index for minimization over the time interval [to,tf] 


r tf 

(4.2) J(u) = ip{y{t f ),t f ) + / L(y(t), u(t),t) dt 

Jt 0 

where the Lagrangian is given by 

(4.3) 


L{y(t),u(t),t) = ~ [y T (t)Qy{t) + u T {t)Ru(t)] . 


The symmetric, nonnegative definite matrix Q and symmetric, positive matrix R weight the state and control 
input while the function penalizes large terminal values of the state. Energy considerations can 

be used to specify both Q and xj). As detailed in [2], an appropriate choice of Q, which arises from the 
minimization of the kinetic and potential energies, is a multiple of the mass matrix. Similarly, the choice 

1 


(4.4) 


V»(y(*/).f/) = ~y (t/)n/y(f/) 


minimizes the final energy when lij is specified as a positive matrix. In the examples which follow, Q and 
11/ were chosen as 


(4.5) 


Q 


d i 


K 


Q 


n, 


^4 


K 


Q 


where K and Q are given in (3.4) and di, • * • , are integer weights. 

The Hamiltonian associated with this system is 

(4.6) H (y, A, u, t) = L(y, u, t) + A T /(y, u, t) 

where A (t) € ]R 2 ( m+1 ) j s the adjoint variable or Lagrange multiplier. It should be noted that the state 
equation (4.1) satisfies 


V = 


8H 

d\ 


where ^ denotes the gradient of H with respect to A. 
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The minimization of (4.2) is constrained by the system (4.1). To pose this as an unconstrained op- 
timization problem, we incorporate the constraints via the Lagrange multiplier and consider the modified 
performance index 


^(w) = ^y T {tf)^fy(tf) + J [L(y,u,t) + \ T {t){f{y,u,t ) -»]] dt 

(4-7) 

= ^y T ( t f)Rfy( t f) + [H(y,u,t) - \ T (t)y] dt. 


The minimum of the constrained functional J occurs at the minimum of the unconstrained functional J 
which in turn occurs when dJ = 0 (see [6, 26, 27]). Enforcement of this condition yields the necessary 
adjoint condition 


(4.8) 


A = 


dH 

dy 


\{t f ) = n f y(t f ) 


and the stationary condition 
(4.9) 



Note that the terminal condition on the adjoint variable is chosen to satisfy the transversality constraint 
for the system. This provides the framework employed in the finite dimensional linear, nonlinear, and 
perturbation control methods discussed next. 


State Tracking Problem 

The goal in many applications entails the control of system dynamics to a specific trajectory s(t) given 
observations 

(4.10) y ob (t) = Cy(t), 

in ]R*, of the state dynamics. An appropriate performance index for this case is 

l*f 

J(u) =if>(Cy(tf) -s(tf),t f ) + / L(y(t) — s(t),u(t), t) dt 

Jt 0 

where L and ip are given by (4.3) and (4.4), respectively. The final time boundary condition for this choice 
is then 


X(t f )=C T U f [Cy(tf) - s(tf)} . 

4.1. Linear Optimal Control. At low drive levels with magnetic biases, experimental data has indi- 
cated a nearly linear relation between input currents to the solenoid and strains output by the transducer. 
This is reflected in the model response and can be employed when designing a control method for such 
regimes. For low drive level applications, reasonable approximate models and control methods can be at- 
tained through linearization about an appropriate input no- One choice is the coercivity value uo = u c at 
which M(u c ) = 0. In this case, the approximate linear control operator B is 

dlid 

Bu(t) = 2M S — — (u c )u(t)b 
ou 
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where 


&M M an M{ rr 

du U kS — OL[M an — Mj, 


+ ncM s 


ici* f— ) + 

a V a / Hl ff 


Under this approximation, the corresponding first-order system is 

y(t) = Ay(t) + Bu(t) + G(t) 

(4.11) 

2/(0) = 2/o • 

The stationary condition (4.9) then yields the optimal control 

u*(t) = -R’ 1 B T \(t) 

while the state constraint (4.11) and adjoint condition (4.8) yields the optimality system 


(4.12) 


y(t) 


. m . 



A -BR~ l B T 
-Q -A T 


y(t) 

m 


+ 


G(t) 

0 


y{t 0 ) = yo 

\{t f ) - u f y(t f ) . 

The construction of the optimal control requires the solution of the two-point boundary value problem (4.12). 
Due to the linearity of the system, however, a fundamental solution matrix can be employed to formulate 
the optimal control as 


(4.13) 


u*{t) = -R~ 1 B 1 p(t)»(i)-r(t)] 


where II (t) solves the differential Riccati equation 


(4.14) 


-n = 4 T n + iL4 - hbr- 1 b t u + q 


n(t/) = n f . 


The perturbation variable r(t) € ]R 2 ( m+1 ) is obtained through integration of the final time system 

r(t) = - [A - Bfi _1 B T n] T r(t) +IIG(i) 
r(t f ) = 0. 

In this manner, the solution of a system with split conditions at the initial and final times is replaced by 
solution of systems with only final time conditions. 

Two special cases are sufficiently common in applications to warrant further discussion. The first con- 
cerns the infinite time problem while the second characterizes Riccati solutions and optimal controls when 
input forces are periodic. It is important to note that in these cases as well as the general finite time for- 
mulation, the control (4.13) acts in a feedback manner on current states of the system. This will not be the 
case with the nonlinear control method. 
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4.1.1. Infinite Time Problems. For the strongly dissipative systems under consideration, it is reason- 
able to assume that {A,B) is stabilizable and (A,C) is detectable (see (4.10) for discussion of the observation 
operator C). In this case, as t — ► oo, y and u approach 0 at a sufficient rate to guarantee the existence of 
the performance index 

(4.15) J(w)=-/ [y T (t)Qy{t) +u T {t)Ru{t)] dt. 

L Jto 

The Riccati matrix used to characterize the feedback control (4.13) is the solution to the steady state 
algebraic Riccati equation 

(4.16) A^ + UA-llBR-^^ + Q = 0 

and is now constant. Similarly, the decay of solutions implies that the perturbation component is given by 
r(t) = 0. Hence implementation of the method requires only the offline solution of a Riccati solution followed 
by online feedback on observed states. 

4.1.2. Periodic Problems. A second case which commonly arises in applications is that in which 
the exogenous force G(t) models periodic or oscillatory inputs to the system. If r denotes the fundamental 
period for all frequencies present, an appropriate performance index is 

J(u) = \ J o [ y T ( t )Qy ( t ) + « T (*)-Ru(t)] dt . 

Under the hypotheses of stabilizability and detectability, it is shown in [3, 12] that the optimal control (4.13) 
can be formulated in terms of the solution to the algebraic Riccati equation (4.16) and the solution to the 
periodic perturbation system 

f(t) = -\a- br~ 1 b t u} t r(t) + ncm 

(4.17) 

r(0) = r(r) . 

This yields a feedback algorithm which is efficient to implement in many applications. 

4.1.3. Numerical Example — No Exogenous Force. To illustrate the performance and limitations 
of the linear control method, we consider the use of the linear control law in the nonlinear system 

y(t) = Ay(t) + [B(u)](«) 

(4.18) 

y(to) = yo 

for various magnitudes of the initial value yo. To obtain these values, the uniform force p(£, x) = go sin(107r£) 
was applied to the uncontrolled beam for to = 0-45 seconds and then terminated. The initial value yo was 
taken to be the state at the time to when control was initiated. Control inputs were computed through 
minimiz ation of the infinite time performance index (4.15) with R = 5 x 10 2 and d\ — d<i = 5 x 10 8 in the 
definition (4.5) for Q. 

The implementation of the method in this manner provides a numerical illustration of effects which may 
be observed when control currents computing using a linear model and control law are fed back into the true 
physical system having nonlinear actuators. While we are not providing here the full convergence analysis 
and model fits to a physical apparatus, numerous experiments have demonstrated the validity of the model 
[9] and the trends illustrated by these numerical results. 
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The uncontrolled and controlled beam displacements at the point x = 3£/5 (see Figure 4) with go = 1 
are plotted in Figure 5b while the displacements generated with go = 100 are plotted in Figure 5d. The 
relationship between the input magnetic field H = ril = nu and magnetization M for the two cases are 
given in Figures 5a and 5d. It is noted that at the low drive level, the relationship between H and M is 
approximately linear and the feedback of the linear control u into the nonlinear system (4.18) is very effective. 
At the higher drive levels, however, the relationship between H and M displays significant hysteresis which 
leads to energy loss and time delays in the input. In this case, the linear control law (4.13) does not provide 
the capacity for accurately quantifying and incorporating the hysteresis and subsequent delays which in turn 
produces the loss in control authority observed in Figure 5d. This illustrates that while the linear control 
method can be effective at low drive levels, it does not provide the accuracy necessary for moderate to high 
drive level applications. For such regimes, control methods which incorporate the actuator nonlinearities are 
required. 



Figure 5. Feedback of linear law (4.13) into the nonlinear system (4.18). Relationship between magnetic 

field H and magnetization M; (a) go — 1 and (c) go = 100. Uncontrolled ( ) and controlled ( ■ »■ ■■■ ) beam 

trajectories at the point x = 3^/5; (b) go = l and (d) go — 100. 
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4.1.4. Numerical Example - Periodic Exogenous Force. A second regime common in structural 
applications is that in which exogenous disturbances are periodic (e.g., oscillating mechanical components). 
In this case, a semidiscretization in the spatial variable yields a system of the form (4.11) where G(t) is 
periodic. To illustrate, the spatial uniform exogenous force g(t,x) = </o sin(107rt) was applied throughout 
the time interval [0, 2.5]. The uncontrolled trajectories at x = 3^/5 with g 0 = 1, 100 are plotted in Figure 6b 
and Figure 6d, respectively. Both cases exhibit a beat phenomenon due to the close proximity of the 5 Hz 
driving frequency with the 6.1 Hz natural frequency for the beam (see Section 2.3). 

Controlling currents were computed via (4.13) with intermediate perturbation solutions obtained through 
integration of the system (4.17). Inputs for the low and high drive level cases are illustrated in Figures 6a 
and 6c. As in the previous example, the linear control method is highly effective at low drive levels where 
the linear model is accurate. At the high drive level in which the actuators are advantageous, however, 
the input exhibits significant hysteresis which acts as a phase delay to the system. The result is a loss in 
control authority to the extent that controlled beam trajectories actually have larger magnitudes than the 
uncontrolled beam. This further motivates consideration of a nonlinear control method. 



(c) (d) 

Figure 6. Feedback of linear law (4.13) into the nonlinear system (4.18). Relationship between magnetic 

field H and magnetization M; (a) p 0 = 1 and (c) g 0 = 100. Uncontrolled ( ) and controlled ( — ) beam 

trajectories at the point x = 3^/5; (b) g 0 = 1 and (d) g 0 — 100. 
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4.2. Nonlinear Optimal Control. We consider here the problem of constructing a nonlinear control 
for the system 


(4.19) 


y(t) = Ay(t) + [B(u)}(t)+G(t) 
y(t 0 ) = yo ■ 


In this case, minimization of the of the performance index (4.2) or (4.7) yields the optimality system 


(4.20) 


y(t) 


Ay(t) + [B(u)](f) + G(t) 

X (t) _ 


_ —A T X(t) - Qy(t) 


y(to) = yo 


K*f) = n /y(*/) 

where the optimal control satisfies 

u-(t) = -R- 1 [BT(u*)}(t)\(t). 

Due to the nonlinear nature of the input operator B(u), decomposition of the system matrices in terms of a 
fundamental matrix solution is not possible which prohibits efficient solution in terms of a Riccati matrix. 
To this end, we consider the approximation of the full two-point boundary value problem (4.20) or the 
equivalent first-order system 

z(t) = F(t,z) 


(4.21) 


E 0 z(t 0 ) = [ 3 / 0 , 0] r 


where z = [y, A] T and 


(4.22) 


Efz(t f ) = [0,11 f y(t f )] T 


F(t,z) 


Ay(t) + [£?(«)](£) + G(t) 
-A T X(t) - Qy(t) 


Eq = 


I 0 
0 0 


E f = 


0 0 
0 I 


Here I denotes a 2(m+ 1) x 2(m-bl) identity matrix where m + 1 denotes the number of basis functions 
employed in the spatial approximation (3.2) of the state variables. 

The solutions to the system (4.21) can be approximated through a variety of methods including finite 
differences and nonlinear multiple shooting. To illustrate a finite difference approach, we consider a dis- 
cretization of the time interval [to, t/] with a uniform mesh having stepsize At and points to, ti, ■ • ■ , tjv = t/. 
The approximate values of 2 at these times are denoted by 20 , • * • , z^. A central difference approximation 
of the temporal derivative then yields the system 

[ z j + 1 — z j] — g z j) + 1 )] 


(4.23) 


EqZq = [y 0 , 0]^ 


for j = 0, • • • , N - 1. 


EfZN = [0,II/2/(f/)] r 
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The determination of a solution vector Zh = [^o, ■ • • , z^\ to (4.23) can then be expressed as the problem 
of finding Zh which solves 

(4.24) ?{zh) = 0. 

For the difference method and boundary conditions considered here, T{zh) £ has the form 


F{ z h) = 


Fo 

Fi 


Fj — ^ [ 2 j+ 1 z j) 2 \F(tj,Zj) 4- F(tj+i, Zj+\)\ 


6 ( zq , zn) = EqZq + EfZjsi 


2/o 


0 

0 


n f y{tf) 


Fn-i 

b{zo,z N ) 

A quasi-Newton iteration of the form z£ +1 = zfc + where ^ solves 
(4-25) F(4)& = -^(4). 

is then used to approximate the solution to the nonlinear system (4.24). The 4(A r +l)(m+l)x4(A r -fl)(m-f-l) 
Jacobian F'(zfc) has the form 

So Ho 

Si R\ 


F\z h ) 


Eq 


Sjv-i Hjv-i 


where 


The matrix «4(£i) is the linearization 


which yields the representation 


Si ~ At 1 2 A ^ 


A( t ,) = d £(k,z<) 


1 

' I 

0 ' 

1 

' A 

&*[<] ‘ 

At 

_ 0 

I 

” 2 

. Q 

-A T 


for Si. The representation for Ri is similar. 

For this application, direct solution of (4.25) is infeasible due to the large number of variables required 
to resolve the solution over a reasonable time interval. The structure of the Jacobian can be employed, 


16 



however, to reduce both memory and computational requirements to the level of solving 4(m+l)x4(m + l) 
systems. To this end, we express the Jacobian in the analytic LU decomposition 

Ftf) = LU 


where 


So 


Si 


L ~ 


S N -i 0 

TV — 2 TV — 1 


E 0 -EoiSo'Rv) ••• EoHi-iyiS-'Ri) E f + E 0 (-lyiS^R,) 


i=l 


i=l 


I Sq 1 Rq 

I S-'Rx 


u = 


I Sn-iRn-i 

I 


The solution of the system (4.25) is then obtained through direct solution of the system lower triangular 
system L(£ = — followed by direct solution of the upper triangular system 


Remark 1: The conditioning of the matrices Si and Ri is partially governed by the choice of state weights 
d\,d 2 in Q (see (4.5)). The conditioning is improved through the choice of values on the order of 10 3 or less. 
To maintain control authority, this dictates the choice of control weights R to be on the order of 10 “ 3 or less. 
For these choices, the component matrices in the lower and upper diagonal systems are well conditioned. 


Remark 2: The spatial approximation of the beam model was fully resolved with m = 12 basis functions 
and the results reported here were obtained with m = 16. This yields a total of 4(m 4- 1)(N + 1) = 13940 
coefficients to be obtained when using a stepsize of At = .01 on the time interval [0.45,2.5]. To test the 
efficiency and memory requirements necessary for extending the method to problems in two space dimensions, 
we also ran the problem with m — 144 basis functions. This would correspond to a discretization with 12 
basis functions in each spatial dimension and yields a total of 118,900 unknowns to be obtained. The 
method is sufficiently efficient so that even in this range, computations could be performed on a workstation 
in Mat lab. 


Remark 3: The matrix products arising in the lower triangular system L can be formed recursively. 
Moreover the component matrices have significant structure which can be utilized when forming the matrix 
products. The utilization of inherent recursions and structure is necessary when considering systems in two 
space dimensions which can have in excess of 500 states. 
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4.2.1. Numerical Example — No Exogenous Force. The use of the nonlinear control method is 
illustrated in the context of the cantilever beam driven for 0.45 seconds by the uniform force g(t, x) = 
100 sin(107rt) at which point the force was terminated and control initiated. The control inputs were com- 
puted using the approximation method (4.23) for the two point boundary value problem (4.21) on the time 
interval [to,tf] = [0.45,2.45]. The control weights were taken to be d\ = = 5 x 10 2 and R = 5 x 10 -4 . 

This yielded an optimal current which was then applied as an open loop control to the system. The resulting 
controlled trajectory at the point x = 3^/5 is compared with the uncontrolled trajectory in Figure 7. The 
corresponding relationship between the input magnetic field and output magnetization is plotted in Figure 8. 
It is noted that the model-based nonlinear control law very adequately incorporates the inherent hysteresis 
in the transducer and provides complete attenuation within 0.5 seconds of being invoked. This illustrates 
the performance of the nonlinear control law and capabilities of the magnetostrictive transducers under ideal 
operating conditions. 

One difficulty with an open loop control law of this type is its lack of robustness with respect to uncer- 
tainties in operating conditions. Such uncertainties can be due to unmodeled dynamics, changing operating 
conditions, or slight delays or phase shifts due to filters, etc., and are present to some extent in all experi- 
mental systems. 

To illustrate the effect of uncertainties on the performance of the open loop control, we consider the 
same system with the control applied 0.03 seconds late. This is a very reasonable scenario in experiments 
and must ultimately be accommodated by the control law. The uncontrolled and controlled trajectories for 
this case are depicted in Figure 9. The slight delay in the initiation of the control input results in a complete 
degradation of control authority (compare with the attenuation in Figure 7 with no delay). This illustrates 
the necessity of feeding back some form of state information and motivates consideration of perturbation 
control methods. 



Figure 7. Uncontrolled and controlled beam trajectories at the point x = 3^/5; (uncontrolled), ■ 

(controlled). 
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Figure 8. Input magnetic field H = nJ and output magnetization M. 



Figure 9. Uncontrolled and controlled beam trajectories at the point x ~ 3£/5 with control initiated 0.03 
seconds late; (uncontrolled), — (controlled). 

4.2.2. Numerical Example — Periodic Exogenous Force. The techniques for computing the open 
loop nonlinear control for systems with exogenous forces are identical to those employed in Section 4.2.1; one 
simply modifies F in (4.21) by the appropriate exogenous force. To illustrate, the force g(t , x) = sin(107r£) was 
applied for the full time interval [0.2.5] with the optimal control computed for the interval [f 0 , tj = [0.45, 2.5]. 
The resulting beam trajectory and inputs are plotted in Figure 10 and Figure 11. A comparison of Figures 
10 and 6 indicates that reductions on the order of those obtained in the low drive level linear case can be 
obtained with the nonlinear law. Figure 11 illustrates that following an initial transient phase, the input 
relation settles into a hysteretic periodic cycle with the frequency matching that of the driving input. 
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In this case, the system is subject to uncertainties in the measured exogenous force in addition to the 
operating uncertainties discussed in Section 4.2.1. This can include perturbations in frequency or phase 
which can destabilize a feedback method and degrade open loop attenuation if unincorporated. In Figure 12, 
we illustrate the trajectories of a beam subjected to the force 

( 100sin(107rt) , t < .45 

<?(<,*) = { j 

( 100sin(147rf — 1.87r) , .45 < t < 2.5 

with the factor of 1.87T included to ensure the continuity of g. The effect of the frequency change can be 
noted in that 12 oscillations are now present in the control interval [.45, 2.5] compared with the 10 oscillations 
noted in Figure 10. The open loop control was computed for the assumed force g{t y x) = 100 sin(107rt) and 
was applied 0.03 seconds late. It is noted that the control attenuation is completely degraded by these 
uncertainties and that further robustness must be incorporated in the method. 

Remark 4: The persistence of beam vibrations in spite of the control input indicates a physical limitation 
of the actuator setup rather than a deficiency in the control formulation. To attain greater attenuation, one 
must investigate controllability issues related to physical criteria such as actuator number and placement. 
The degree to which such physical issues play a role depends upon the application and in many cases, 
attenuation on the order of that observed in Figure 10 is sufficient. 



Figure 10. Uncontrolled and controlled beam trajectories at the point x = 3^/5; (uncont rolled), — 

(controlled). 
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0 0.5 1 1.5 2 2.5 

Time 

Figure 12. Uncontrolled and controlled beam trajectories at the point x = 3^/5 with control initiated 0.03 
seconds late and a 2 Hz frequency perturbation; (uncontrolled), 1 (controlled). 

4.3. Perturbation Control. As illustrated in the last example, a purely open loop control law suffers 
from lack of robustness with regard to system uncertainties. For various classes of uncertainties, robustness 
can be significantly enhanced through consideration of perturbation control techniques [6, 27]. In these 
methods, the system is linearized about the optimal control pair (u*(£), y*(t)) obtained through solution 
of the two-point boundary value problem (4.20) or (4.21). A feedback control <5u*(£) is then designed to 
attenuate perturbations in the system due to uncertainties in the exogenous force or uncertainties in initial 
conditions as depicted in Figure 13. Both are common in applications with perturbed initial conditions often 
due to uncertainties in the starting time for the open loop control. Because LQR theory can be invoked to 
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Figure 13. Optimal open loop controlled state and neighboring state due to perturbed initial condition. 


construct the perturbation control the implementation of the method is very efficient once the open 

loop control pair has been determined. 

The perturbation control system can be obtained by expanding the augmented cost functional and 
constraint equations through higher-order terms and employing the simplifications provided by the fact that 
u*{t) and y*(t) minimize the first-order optimality system. Since d J = 0 for the optimal pair (ti*(£), y*(t)), 
expansion of the augmented cost criterion (4.7) through second-order terms and constraints to first-order 
yields 


(4.26) 


and 



Hyu 

H uu 


5y 

Su 


dt 


(4.27) 


Sy(t) = ASy(t) + B6u(t) + SG(t) 
Sy{ 0 ) = yo 


where 5u and Sy are first-order variations about u* and y*. The optimal perturbation control 6u* is that 
which minimizes (4.26) subject to (4.27). 

For the Hamiltonian (4.6), the second variation S 2 J is given by 


(4.28) 


S 2 J = 


5 I'' {(QSy,Sy) + (RSu,Su)} 
1 Jto 


dt 


so that the LQR theory outlined in Section 4.1 can be directly employed to obtain and 5y*(t). The 

overall control for the system is then taken to be u*(t)+6u*t(t) with the optimal state given by + 

For implementation purposes, it should be noted that the optimal open loop control u*(t) can be computed 
offline leaving only Su*(t) to be computed online. 


4.3.1. Numerical Example — No Exogenous Force. The performance of the method is first illus- 
trated in the context of Example 4.2.1 in which a perturbed initial value is introduced through application 
of the optimal control u*(t) to the system 0.03 seconds late. As noted in Figure 9, this perturbation is 
sufficient to destroy the authority of the open loop control. To accommodate these perturbations, we employ 
the control law 


5u*(t) = —R~ 1 B T USy(t) 
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where II satisfies the algebraic Riccati equation (4.16). The resulting control trajectory is illustrated in 
Figure 14. It is noted that through the use of the feedback perturbation control, attenuation comparable to 
that for the unperturbed system (see Figure 7) is obtained. This provides a significant enhancement of the 
method with respect to perturbations in initial conditions. 



Figure 14. Uncontrolled and controlled beam trajectories at the point x = 3t/5 with control initiated 0.03 
seconds late; (uncontrolled), — (controlled). 


4.3.2. Numerical Example - Periodic Exogenous Force. Systems driven by an exogenous force 
are subject to force perturbations in addition to initial uncertainties or delays in control implementation. As 
illustrated in Example 4.2.2, perturbations from the expected 5 Hz force to a measured 7 Hz force completely 
degrade the open loop control. In this case, 

Sg(t,x) — 100[sin(147rt — 1.87t) — sin(107rt)] 
over the time interval [0.45, 2.5], and the perturbation control has the form 

5u*(t) = —R~ l B T [ 7 r6y*(t) — <5r(£)] 


where Sr(t) solves 

Sf(t) = - [A - BR- 1 B t U\ T 6r(t) + II 6G(t) 

6r(0) = Sr(r) . 

In addition to the force perturbation, a perturbed initial condition due to a 0.03 second delay in control 
initiation was included in the system. 

The uncontrolled and controlled beam trajectories at the point x — 3£/5 are compared in Figure 15. 
It is noted that while the trajectories differ in frequency due to the combined open and closed loop effects, 
significant attenuation is attained throughout the time interval due to the feedback perturbation control 
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component. A comparison with Figure 10 indicates that the controlled trajectory is comparable in magni- 
tude to that in the perturbed case even though the uncontrolled displacement is significantly larger. Such 
attenuation levels have also been noted with larger frequency perturbations (e.g., a perturbed driving fre- 
quency of 18 Hz). Hence the perturbation control provides a feedback methodology which is highly robust 
as well as efficient to implement. 



Figure 15. Uncontrolled and controlled beam trajectories at the point x = 3^/5 with a 2 Hz force pertur- 
bation and control initiated 0.03 seconds late; (uncontrolled), — — (controlled). 

5. Concluding Remarks. This paper addressed the development of a physics-based control method- 
ology appropriate for magnetostrictive actuators in moderate to high drive level regimes. At such drive 
levels, these materials exhibit significant hysteresis and nonlinear dynamics which must be incorporated in 
the model and control method to attain the full potential of the actuator (both experiments and numer- 
ical simulations have demonstrated that linear methods fail at such drive levels). For various structural 
applications, it is also necessary to control both transient and steady state dynamics. 

To attain these objectives, a model based upon ferromagnetic mean field theory was used to characterize 
the actuator dynamics including the inherent hysteresis and nonlinearities. This provided a method of accu- 
rately quantifying multiple frequencies and transient dynamics. Optimal control theory was then employed 
to obtain an open loop control which incorporated the actuator hysteresis and nonlinearities. This nonlinear 
control was combined with a perturbation feedback control to attain a hybrid method which was highly 
robust and efficient to implement. Finally, the efficacy of the method was demonstrated through numerical 
examples. 

We note that the method described here does not address the minimum time control problem nor 
does it actively enforce admissibility criteria. If time minimization is desired, the control problem can 
be reformulated with the final time and final adjoint values treated as components of the solution. For 
applications which require that the control u(t) fie in an admissible region, the stationary conditions (4.9) 
must be replaced by some form of the Pontryagin minimum principle. Details concerning both cases can be 
found in [6, 27]. 
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In its present form, the method is currently designed for linear structural models. While it was illustrated 
in the context of a PDE-based thin beam model, the flexibility for employing large discretization limits in 
Matlab (in excess of 144 basis functions) indicates that the method can be directly applied to certain linear 
plate and shell models. For larger problems in which the number of spatial variables or time steps prohibits 
global optimization over the full time interval, piecewise methods of the type described in [20, 21] can be 
employed to obtain suboptimal solutions over each time step. These piecewise states and controls can then be 
patched together to obtain a global solution over the full time interval. Finally, we note that the extension of 
this method to nonlinear structural models is also important due to the advantages of high output actuators 
in such regimes and is under current investigation. 
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